Kinetics of Phase Separation in Thin Films: Lattice versus Continuum Models for Solid 

Binary Mixtures 
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A description of phase separation kinetics for solid binary (A,B) mixtures in thin film geometry 
based on the Kawasaki spin-exchange kinetic Ising model is presented in a discrete lattice molecular 
field formulation. It is shown that the model describes the interplay of wetting layer formation and 
lateral phase separation, which leads to a characteristic domain size £(t) in the directions parallel 
to the confining walls that grows according to the Lifshitz-Slyozov i 1 / 3 law with time t after the 
quench. Near the critical point of the model, the description is shown to be equivalent to the standard 
treatments based on Ginzburg-Landau models. Unlike the latter, the present treatment is reliable 
also at temperatures far below criticality, where the correlation length in the bulk is only of the order 
of a lattice spacing, and steep concentration variations may occur near the walls, invalidating the 
gradient square approximation. A further merit is that the relation to the interaction parameters 
in the bulk and at the walls is always transparent, and the correct free energy at low temperatures 
is consistent with the time evolution by construction. 



I. INTRODUCTION 

A basic problem of both materials science [JE! and sta- 
tistical mechanics of systems out of equilibrium dHHH is 
the process of spinodal decomposition of binary (A,B) mix- 
tures. When one brings the system by a sudden change of 
external control parameters (e.g., a temperature quench) 
from an equilibrium state in the one-phase region of the 
mixture to a state inside of the miscibility gap, thermal 
equilibrium requires coexistence of macroscopically large 
regions of A-rich and B-rich phases. Related phenomena 
also occur in systems undergoing order-disorder phase tran- 
sitions, and it is also possible that an interplay of ordering 
and phase separation occurs, particularly in metallic alloys 
with complex phase diagrams P, 0, Q . 

The recent interest in nanostructured materials and thin 
films has led to consider the effects of walls or free sur- 
faces on the kinetics of such phase changes 0, H, H, El EH, 
El, El EI El El El El El- In a binary mixture, it is 
natural to expect that one of the components (say, A) will 
be preferentially attracted to the surface. Already in the 
one phase region, this attraction will lead to the forma- 
tion of surface enrichment layers of the preferred compo- 
nent at the walls, but the thickness of these layers will be 
small, viz., of the order of the correlation length £ of con- 
centration fluctuations in the mixture [2(|. Only near the 
critical point this surface enrichment becomes long ranged 
[UHl]. For temperatures below the critical point, how- 
ever, the behavior is more complicated due to the inter- 
play between bulk phase separation and wetting phenom- 
ena [H, [H [H, [H, [23, [H HI ■ In a semi-infinite geometry, 
a B-rich domain of a phase separated mixture would be 
coated with an A-rich wetting layer at the surface, if the 
temperature is above the wetting transition. Of course, for 



thin films the finite film thickness D also constrains the 
growth of wetting layers: E.g., for short range forces be- 
tween the walls and the A-atoms the equilibrium thickness 
of a "wetting layer" is of the order of £ln(D/£) [HHHH, 
while at lower temperatures where the surface is nonwet the 
thickness of the surface enrichment layer again only is of 
the order of £. Also the phase diagram of the system in a 
thin film geometry differs from that of the bulk, in analogy 
to capillary condensation of fluids [H [3l|, [H, [H, HI] , and 
very rich phase diagrams may occur, in particular, if the 
surfaces confining the thin film are not of the same type 
[H [H, [H H3] ■ These changes in the phase behavior of 
thin films are reflected in the kinetics of phase separation 
in these systems 0, El, of course. 

Despite the large theoretical activity on surface effects on 
hase separation kinetics @, ©, El El El, El, El El El, 
the applicability of these works on experiments is very 
restricted. Actually, most of the experiments deal with 
thin films of fluid binary mixtures \7 ,11,16] 
all theoretical works [1 M El, El, [3 El, (S , _ 
with "model B" [11] , where hydrodynamic interactions are 
neglected, and hence this model is really appropriate only 
for solid mixtures. 

A second restriction on the applicability of the theory is 
the fact that it is based on the time-dependent Ginzburg- 
Landau model p], H, H 0, H, H, and thus is applicable 
only when the correlation length £ is very large, i.e. in 
the immediate vicinity of the bulk critical point. For a 
thin film with short range surface forces, the "standard 
model" is based on a free energy functional of an order 
parameter ip, and this functional consists of a bulk term 
(Ft,) and two terms representing the two surfaces 51, 52 
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Here ?/>(r) is the order parameter which is proportional to 
the density difference between the two species; it is normal- 
ized such that the coexisting A-rich and B-rich bulk phases 
for temperature T < T c b (X^ being the bulk critical tem- 
perature) correspond to if) = ±1, respectively. All lengths 
are measured in units of 2£, with £ denoting the bulk corre- 
lation length at the coexistence curve. The terms Fsi and 
Fs2 are the local contributions from the surfaces 51 and 
52, which for a film of thickness D, are located at z = 
and z = D, respectively, orienting the z-axis perpendicu- 
lar to the surfaces, while p denotes the (d — 1) coordinates 
parallel to the surfaces, d being the dimensionality. In Fsi 
there are parameters g, 7, and hsi, which can be related to 
the temperature T and various parameters of an Ising fer- 
romagnet with a free surface at which a surface magnetic 
field Hsi acts S, HI 
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In Eqs. (HI [5]), a simple cubic lattice was assumed, the z- 
axis coinciding with a lattice axis. Nearest neighbor Ising 
spins in the lattice interact with an exchange coupling J, 
except for the surface plane at z = where the exchange 
coupling is J5. For deriving Eqs. ([1] [2]) from a layerwise 
molecular field approximation, the limit £ — > 00 needs to 
be taken, with ip(p, z) — m n (p, t)/mb, where n is the layer 
index of the lattice, n — 1 being the surface plane, and nib 
[= V3(l - T/T ch ) 1/2 } the bulk magnetization (note that 
Wb — > in the considered limit). Finally, Fs2 describes 
analogously, the surface free energy of the surface at z = D, 
with h$2 being related to the surface field Hs2 analogously 
to Eq. In Eq. ©, £ is measured in units of the lattice 
spacing a of the molecular field lattice model, which often, 
for the sake of convenience, we set to 1. 



Dynamics is associated to the model assuming that in the 
bulk the order parameter evolves according to the (nonlin- 
ear) Cahn-Hilliard equation p], 0, H, 0, IE [y| 



d - - 

—^(r,T)=-V-J(r,T) = V 
or 
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(6) 



while the surfaces amount to two boundary conditions each 
20], which can be written as 



T -i>(p,0,T) = -^r^y 



= h sl +gip(p,0,T) + j—\ z=o , (7) 
oz 



Jz(p,0,r) 



d_ 

dz 



0. (8) 



Equation ([7]) describes a non-conserved relaxation ( "model 
A" (38|) for the order parameter at the surface, To setting 
the time scale. Since ip(p,0,T) relaxes much faster than 
the time scales of phase separation away from the surface, 
one may put To = 0. The equations for z — D are fully 
analogous to those for z — 0. We emphasize that Eq. ([6]) 
can also be derived (4fj| from a continuum approximation 
to a molecular field approximation to a description of a 
Kawasaki kinetic Ising model [4l|, and Eqs. Q [5]) can be 
derived from a Kawasaki model with a free surface as well 
[2fJ. However, it is clear that the model, Eqs. ©El [5]), can 
only represent the molecular field Kawasaki kinetic Ising 
model accurately when £ 3> one lattice unit. This fact was 
already discussed in [20| in the one-phase region, where 
the linearized molecular field equations on the lattice were 
solved to describe the kinetics of surface enrichment for 
T > T c b, and it was found that the lattice and continuum 
theories agree when approximations such as exp(— 1/£) — 
1 — l/£ become valid. 

Thus, the model Eqs. ©El [5]) can describe a solid binary 
mixture accurately near the bulk critical point only: far 
below T c b, terms of order if) 6 and higher would be needed 
in Eq. |T]) already, and when £ is of the order of the lattice 
spacing a, also higher order gradient terms (V 2 ip) 2 etc. 
would be required for an accurate continuum description. 
However, the numerical solutions of Eqs. (JH1 [3 [B]) require 
anyway a discretization: Often a mesh size Ax = Ay = 
Az = 1 is used to solve these equations @, [E IE IE 
IE IE IE] • This essentially corresponds to a lattice with 
spacing 2£, rather than the spacing a of the underlying 
Ising lattice. Having in mind applications of the theory 
at temperatures that are not close to the critical point, 
however, £ is of the order of 1. 

Thus, it is plausible that one could solve with a compara- 
ble effort the original (nonlinear) molecular field equations 
for the Kawasaki kinetic Ising model on the lattice, that 
are underlying this theory. The advantage of such a treat- 
ment clearly would be that the model has a well-defined 
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microscopic meaning at all temperatures. Near the criti- 
cal temperature, the results of this treatment should be- 
come indistinguishable from the solution of Eqs. ([HJ[71[8]), 
of course. 

The purpose of the present work is to show that indeed 
such a lattice mean field approach to spinodal decompo- 
sition in thin films (and in the bulk) is both feasible and 
efficient. In Sec. II, we present the discrete lattice analogue 
of Eqs. ||TJ |3J), and (E3 HJ), following up on the work 
by Binder and Frisch [20]. In Sec. Ill, we present various 
numerical results for deep quenches (i.e., temperatures far 
below criticality) and discuss the resulting structure forma- 
tion for a symmetric film, both in directions parallel and 
perpendicular to the walls. In Sec. IV, we present a com- 
parison to the continuum approach, Eqs. ©El [5]). Finally, 
Sec. V summarizes our paper, containing also an outlook 
on simulations where the basic atomistic model is not a 
lattice model, as is the case for spinodal decomposition in 
fluid binary mixtures. 



II. MOLECULAR FIELD THEORY FOR THE 
KAWASAKI KINETIC ISING MODEL IN A THIN 

FILM 

The system that is considered is the ferromagnetic Ising 
model with nearest neighbor exchange on the simple cubic 
lattice, described by the Hamiltonian 

H = —J SiSj — J s SiSj 

<i,J> <i,3> 
bulk 51,52 

-H J2 Si - Hsi S i - H S2 ]T S t . (9) 
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Here Si = ±1, lattice sites are labeled by the index i, and 
the first sum runs over all nearest neighbor pairs except 
those in the surfaces 51 and 52. Note that now for a film 
of thickness D, the surfaces 51 and 52 are located at n = 1 
and n = n m ax = D + 1, n being an index labeling the lat- 
tice planes in the z-direction perpendicular to the surfaces. 
Also note that all distances will be measured in units of 
the lattice constant a. The term —Hj^Si describes the 

i 

Zeeman energy, H being the bulk magnetic field, and the 
sum runs over all sites of the lattice. It should be remem- 
bered that in the interpretation of the Ising model related 
to binary mixtures, spin up corresponding to A and spin 
down to B, H would correspond to a chemical potential 
difference between the species. 

The molecular field equations for the local magnetization 
m n (p) = (Si) (we denote the index i of a lattice site by the 
index n of the plane to which it belongs and a coordinate 
p in this plane) become [20| 



m n (p) = tanh 



J 



y{p) 



y{p) 



^2m n {p + Ap) 



H 



2 < n < D, (10) 



A/7 



TOi(p) = tanh 



J 



H sl +H 



n 2 {p) + -f ^mi(p + Ap) 

Ap 



, J, n=l, (11) 
m D +x{p) = tanh — — mp(p) + m D+ i(p + Ap) 



k B T 
Hs2 + H 



J 



J 

Ap 

n = D + 1. 



(12) 



Here Ap is a vector connecting site i with one of its 4 
nearest neighbors in a layer. In the bulk, Eq. (Ti"Q|) reduces 
to the well-known transcendental equation 



m b = tanh[-4-(6m b + H/J)], 
kb-1 



(13) 



from which T c b = 6 J/k B straightforwardly follows. 

We also note that Eqs. (|10I11I12[) correspond to the free 
energy 



D+l 



with 

S n 



F=J2(En~ TS n 



1 + m n (p) \ ( 1 + m n (p) 



(14) 



1 - m n (p) ^ (1 - m n {p) 



)}• 



(15) 



E n = - \m n (p)H + — m n (p)(mn_i(p) + m n +i(p)) 
p 

+-Jm n (p)J2 m n(p+ A P)}, 2<n<D, (16) 

Ap 

E x = -Y,{rni(p){H + H x ) + ^j mi (p)m 2 (p) 
p 

+ ij s m 1 (/?)^m 1 (p + Ap)}, n=l, (17) 

Ap 

Ed+i = - {m D+ i(p)(H + H S 2) + ^m D+1 (p)m D (p) 
p 

+^m D+1 (p)^(p + Ap)}, n = D+l. (18) 

Ap 

Of course, Eqs. (TTUl EB E2) can be derived from Eq. (H]) 
via 



\dm„(p) ) 



dm n (p) J T,H,H u {m n (p)y 



= 



(19) 
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The prime on {m„(/j)}' indicates that all local magneti- 
zations are held constant except the one that appears in 
the considered derivative. Of course, in equilibrium there 
is no explicit dependence on p, although there clearly is a 
dependence on n. 

However, already in equilibrium the p-dependence is use- 
ful, when it is understood that the bulk field has a p- 
dependence. Considering a wave vector dependent bulk 
field, one then can derive from Eq. pT7|) for D — > oo the 
wave vector dependent susceptibility By linear re- 

sponse to the wave vector dependent field one finds, as is 
well-known, 



X(k) = 



1 



I -ml 



k B T l-[J(k)/k B T\(l-m*y 



(20) 



where J(k) is the Fourier transform of the exchange inter- 
action, i.e., in our case 



J(k) = 2 J[cos(k x a) + cos(k v a) + cos(fc z a)] ~ 6 J 
From Eqs. (|20l [2~Tj) we readily see that 

r knT 



Jk 2 a 2 . 
(21) 



x(fc) = x(0)/[i + ^V 

and 



x(o) = 1/ 



J 



1 



- k B T r] 



1 — mi 



knT r i 



(22) 
(23) 



Equation (j2"3"|) reduces to the result quoted in Eq. §5§ when 
T — > T c b, but also shows that £ — > when T — > 0. 

As discussed in [20|, we associate dynamics to the Ising 
model via the Kawasaki spin exchange model [4l[ , consid- 
ering exchanges between nearest neighbors only. Follow- 
ing the method of |4(J, one obtains with the help of the 
Glauber [52] transition probability a set of coupled kinetic 
equations for the local time-dependent mean magnetiza- 
tions (Si(t)) = m n (p,t) from the (exact) master equation 
PH in molecular field approximation. Denoting the time 
scale in the transition probability as 7~s, this set of equa- 
tions is: 

(i) 3 < n < D - 1 (bulk case) 



2ts-t: m n {p,t) 
at 



-6m n (p,t) 
+m n -i(p,t) + m n+1 (p,t) 
+ y^m»(p + Ap,t) 



Ap 



[1 - m n (p,t)m n -i(p, t)]tanh- 



J 



k B T l 

i-i(p, t) +^ m n (p + Ap, t) 



f-n+l 



(P,t) 



-m n {p,t) - m„_ 2 (p,t) - 7^ TO n _i (p + Ap, t)] 

Ap 

+ [1 - m n (p, t)m n+1 (p, t)] tanh — — [m„ + i(p, t) 

k B l 

+m n -i(p, t) + ^ m n(P + Ap, t) - m n (p, t) 

Ap 

-m n+2 (p,t) - y^m n+ i( / o + Ap,t)] 

Ap 

+ 2^[1 - m n (p, t)m n (p + Ap, t)} 

Ap 

tanh — —[m n +i{p,t) +m n - 1 (p,t) + 
KB -l 

2J m n (p + Ap', t) - m n +i(p+ Ap,t) 

Ap' 

-m n -x(p + Ap,t) 

-^m„(p + Ap + Ap',t)}. (24) 

Ap' 

Factors such as [1 — m n m n _i] arise from the mean 
field approximations to factors (1 — SiSj) that ex- 
press the fact that an exchange of Sj with Sj changes 
the magnetization Sj only if Sj and Sj are oppositely 
oriented. In Eq. (|2~4")1 , exchanges of a spin at site p 
in layer n with spins in layers n — 1, n + 1, and the 
same layer n need to be considered. In the argument 
of the tanh functions, the difference of the effective 
fields acting on the spins that are exchanged is found. 
One can verify that Eq. (|24|) reduces to Eq. (fTO]) if 
dm n (p,t)/dt = is assumed: in equilibrium, the ex- 
change of a spin with another one is exactly compen- 
sated by the inverse process. The kinetic equations 
near the wall are similar; one has to consider that 
in layer 1 a field H$i is acting, and that no spin ex- 
change into the layer n = (the wall) is possible. 
Hence, for 



(ii) n = 2 



2t s — m 2 (/9, t) = -6m 2 (p, t) + m x (p, t) 
dt 



+7713 (p, i) + ^ m 2 (p + Ap, t) 

Ap 



J 



+[l-m2 {p, t)mx (p, t)} tanh — — [m 3 (p, t) 

k B l 

+7711 (p, t) + ^ m 2 (p + Ap, t) 



Ap 



Hsx 



- -m 2 (p,t) - -y 5Ztoi (p + Ap,t)] 



J 



Ap 



J 



Ap 



+[1 777/2 (p, t)777 3 (p, t)] tanh — — [777.3 (p, t) + mi (p, i) 

Kb J 
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-^m 2 (p + Ap,t) 



Ap 

-m±(p,t) - m 2 (p,t) - y^m 3 (p + Ap,t)] 

Ap 

- m 2 {p,t)m 2 {p + Ap,*)] 

Ap 



tanh — — [m 3 (p,t) + m 1 (p,t) + V m 2 (p + Ap, t)] 

Ap 

-m 3 (p + Ap,t) -m 1 {p+ Ap. t) 

-]Tm 2 (p + Ap + Ap',t), (25) 
Ap' 

and for 



(iii) n = 1 (now only 5 neighbors are available for an 
exchange) 



2T S —mi(p,t) = -5mi(p,t) + m 2 (p,t) 
at 



Y,mi{p + Ap,t) 

Ap 



J 



[1 - mi (p, t)m 2 (p, t)] tanh — — [m 2 (p, t) 

Kb J 



J J 



^mi(p + Ap,i) 



Ap 



-m 3 



(p, t) - mi (p, i) - ^2 m 2 (P + A P> *)] 



^[l-mi(p,t)mi(p + Ap,t)] 



Ap 



tanh [m 2 (p, i) + -y ^ mi (p + Ap ', t) 



knT 



Ap- 



~m 2 



(p + Ap,t) - -£^ mi (p + Ap" ' + Ap,t)]. (26) 



Ap-' 



The equations for n = D and n = D + 1 are analogous. 
The time t in Eqs. (HU [25l [26]) is related to time r in 
Eqs. ©El ED via 



(27) 



T ch /T - 1 



The numerical solutions of the set of equations 
for quenching experiments from infinite temperature to the 
states inside the miscibility gap is the subject of the next 
section. 
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FIG. 1: (a) Layerwise order parameter ^ a v{n) plotted vs. layer 
index n for four different times for the choice D = 29, L = 
128, and Hsi = Hs2 = 1-0. The continuous lines are cubic 
interpolations to the original data, used as guides to the eye. 
(b) Cross-sectional snapshot pictures of the same systems as in 
(a), displaying the magnetization configuration in the xz-plane 
for y = L/2. If the magnetization at a lattice site is positive, a 
black dot is printed, (c) Same as (b), but for a plane parallel 
to the walls, at n — 15. 
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FIG. 2: Layerwise average order parameter ^^(n) plotted vs. 
n for four different times for the cases (a) D = 19, L = 128 and 
(b) D = 9, L = 128. Always Hsi = H S2 = 0.1 is chosen. 
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III. NUMERICAL RESULTS FOR PHASE 
SEPARATION FOLLOWING DEEP QUENCHES 



FIG. 3: Cross-sectional snapshot pictures of the system with 
D = 9, L = 128, H S1 = H S 2 = 0.1 for (a) the xz-plane and (b) 
the zy-plane at n — 5 for the same system as in Fig. [2jb). 



In this section we shall present selected results from the 
numerical solutions of Eqs. ([Ml [23 l^rj]) . choosing three 
values of the thickness of the film (D = 9, 19, and 29, 
corresponding to n m ax = 10, 20, and 30 lattice planes, 
respectively), and a quench at time t = from infi- 
nite temperature to T/T c b = 2/3, for the special cases 
Hsi = H$2 = 1-0, 0.1, and Js = J- The initial conditions 
for m n {p,t) is chosen by taking m„(p, 0) from a random 
uniform distribution between —1 and +1, with the total 
magnetization in the thin film zero. 

At first sight a quench to T/T c b = 2/3 does not look 
like a particularly deep quench. However, one must keep 
in mind that in order to have £ larger (or equal) than a lat- 
tice spacing one must have T > 5.57J/fcB, i.e. much closer 



to T c b = 6J/fes. At the present temperature k^T/J = 4.0, 
the correlation length is as small as £ ~ 0.33 lattice spac- 
ings. Following the standard reasoning for the simulation 
of the Ginzburg-Landau model, Eqs. ([HI L3 0), one should 
choose 2£ as the size of the spatial discretization mesh: 
dealing with the above physical film thickness would re- 
quire rather huge lattices. In addition, Eqs. HI I2I3[) do 
not represent the actual free energies of the lattice model 
[Eqs. (HI EH [H El QUI] accurately at such low tempera- 
tures either. 

Choosing the time unit t s = 1 in Eqs. ([2~4l [231 121))) we 
have found that accurate numerical solutions of these equa- 
tions result already when one chooses a rather large discrete 
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FIG. 4: (a) Layerwise correlation function C n (p,t) plotted vs. 
p, for the choice H S i = H S 2 = 0.1, D = 9,L= 128, t = 10000. 
Data for n = 1,2,3,4 and 5 superimpose almost exactly, (b) 
Fourier transform S n (k,t) of C n (p,t), again resolved with re- 
spect to individual layers. In (b) the dashed line corresponds 
to the Porod tail fc~ 3 . 



time step, 5t = 0.1. Other than this discretization of time 
no approximations whatsoever enter the numerical solu- 
tion. Of course, one always has to deal with a finite sys- 
tem geometry also in lateral directions: unless otherwise 
mentioned we choose L x = L y = L = 128 for all values 
of D and apply periodic boundary conditions in both the 
directions. 

Figure [1] shows data for a typical time evolution for D = 
29 and Hsi = H$ 2 = 1.0. In Fig. Ufa), we show the lay- 
erwise average order parameter, $ av (n) = L~ 2 m n(p, t), 

p 

as a function of the layer index n. One sees that in the 
surface planes the magnetization takes its saturation value 
rather fast, as expected due to the large surface fields. 
Then the magnetization decreases very rapidly, already for 



10 1 10 2 t 10 3 10 4 

FIG. 5: (a) Same as Fig. |4{a) but only for n = 5, and four 
times as indicated, (b) Characteristic domain length £ n (i) of 
individual layers, shown on a log-log plot versus time t. The 
dotted line corresponds to LS growth law t ly/3 . 

n = 3 the magnetization for t — 50 is strongly negative. 
For t — 50 the curve ^> av (n) then exhibits the oscilla- 
tions typical for "surface-directed spinodal decomposition" 

s s m m m m, o, m mm, m in, with a sec- 
ond maximum at n = 7 and a (weak) third one at n = 13. 
Because of the symmetric surface fields, it is expected that 
the profile should be symmetric around the center of the 
film, 

*av(«)=*av(i5-n + 2). (28) 

But in reality this is not obeyed because of finite system 
size and lack of averaging over sufficiently large number 
of independent initial configurations (in our case averaging 
was done only over 5 independent random initial config- 
urations). However, all plots for ^ av (n) have been sym- 
metrized by hand by taking advantage of property (f2"g)) . 
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One can further see from Fig. [TJa) that the thickness 
of the surface enrichment layers at the walls slowly grows 
with increasing time, which is also obvious from Fig. [T^b) 
where we present the snapshot pictures of vertical sections 
(xz plane) . At a later time the position of the second peak 
of ^> av (n) has moved towards the center (it now occurs 
at n = 11 for t = 500) with a pronounced minimum in 
the center of the film. However, for t > 2000 this second 
peak has merged with its mirror image, i.e., now in the 
film center there is a maximum of 4' av (n) rather than a 
minimum. At very late time, this central maximum has 
disappeared again (t = 10000) and the profile looks like 
that of a simple stratified structure. Actually this is not 
the case, as seen in Fig. [ljc), where we have shown the 
snapshot pictures parallel to the surfaces (xy plane) for 
n = 15. As is also observed in the Ginzburg-Landau studies 
of surface-directed spinodal decomposition 0, [jj II 
H3L [li . [I7L HH, in the lateral directions one can observe 
initially a rather random pattern which rapidly coarsens 
with increasing time. 

Figure [5] shows other examples where we have chosen 
thinner films (viz., D — 9 and D — 19) and a much weaker 
boundary field (Hgi = Hs2 = 0.1). Now the amplitude 
of the variation of ^ av {ri) is much smaller, and at late 
times the order parameter profiles across the film have al- 
most no structure. The explanation for this behavior is 
seen in Fig. [3] where we show snapshot pictures of the 
states evolving for the choice D = 9 of Fig. [2Kb): The 
system develops towards a two-dimensional arrangement 
of columns of positive magnetization connecting the two 
walls, and thus in each plane (n — const), there is only a 
weak excess of magnetization in any one direction. These 
results qualitatively do not differ from the numerical stud- 
ies of spinodal decomposition in thin films based on the 
Ginzburg-Landau equation, such as [lH; however, the ad- 
vantage of the present treatment is that the parameters of 
the model have an immediate and straightforward phys- 
ical meaning. The fact that in the late stages there is 
almost no nontrivial structure across the film is also ev- 
ident from a comparison of the pair correlation function 
C n (p,t) [= (m n (Q,t)m n (p,t)) - {mn(d,t)){m n (p,t))] and 
their Fourier transforms S n (k 7 t) in different layers, shown 
in Fig. HI 

In Fig. EJa), the plot of the time evolution of C n (p,t), 
for n = 5, for the same system as in Fig. |3l clearly re- 
flects the coarsening behavior. Note that the apparent 
nonscaling behavior of C n (p,t) is due to strong fluctu- 
ation of the layerwise average order parameter at early 
time. In Fig. [Hb) we plot the layerwise average domain 
size £ n (t) as a function of t extracted from the condition 
C n (p = £ n (t),t) — C„(0,t)/2. The characteristic length 
£ n (t) initially grows rather slowly, for times 10 < t < 1000 
there is considerable curvature on the log-log plot, while 
for t > 1000 the behavior is already close to the standard 
Lifshitz-Slyozov (LS) 0] £(t) oc t 1 / 3 law. The fact that the 
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FIG. 6: Plot of finite-size critical temperature T C (D) versus film 
thickness D, (a) on a linear scale and (b) on a log-log plot to 
display the asymptotic behavior T c b — T C (D) oc D~ 2 of mean 
field theory. The dashed line in (a) marks the bulk critical 
temperature, whereas the continuous line in (b) corresponds to 
the asymptotic theoretical prediction. 



"effective exponent" d[ln£(£)]/<i(ln£) approaches 1/3 from 
below is quite reminiscent of Monte Carlo simulations of 
coarse ning in the two-dimensional Kawasaki spin-exchange 
model [44j, of course. It should be noted that Monte Carlo 
simulations include thermal fluctuations that are absent 
in our molecular field treatment. However, it is generally 
believed [J] that thermal fluctuations are irrelevant dur- 
ing the late stages of coarsening. In view of these facts, 
the similarity of our results with the previous Monte Carlo 
studies of coarsening on lattice models is not unexpected. 
The advantage of the present approach in comparison with 
Monte Carlo, however, is that much less numerical effort is 
needed. 
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FIG. 7: Comparison of the layerwise average order parameter 
profiles across the film for three times (t = 50, 500 and 10000, 
respectively), according to the lattice model (open symbols) 
and the GL model (filled symbols), for (a) D = 9, and (b) 
D = 19, respectively. All data refer to the choice Hsi = Hs2 = 
0.1, k B T/J = 5.875. 



IV. THE CRITICAL REGION: A COMPARISON 
BETWEEN THE LATTICE APPROACH AND THE 
GINZBURG-LANDAU TREATMENT 

Near the critical point one must take into account that 
the critical point is slightly shifted to lower temperatures 
for films of finite thickness as compared to the bulk [2l|, [3l|, 
HI, S| • In Fig. [6^a), we plot the critical temperatures, 
T C (D), for films of finite width as a function of width D. 
The numerical values for T C (D) were obtained by solving 
Eqs. (UHl QU [12]) for H = H SI = H S2 = starting with 
the assignment of uniform magnetization to the layers in a 
random fashion so that the total film magnetization is zero. 
Note that as D — > oo, T C (D) is expected to approach its 



3-ei mean field value 6J//cb = T c b and in the limit D — > 0, 
we expect the 2-d value 4J/fce. In Fig. [UJb), we present 
the deviation of T C (D) from T c h as a function of D on a 
log- log plot. From finite-size scaling theory, this difference 
should vanish as T cb - T C (D) ~ D~ 2 . 

While for D = 29 we have T C (D) ~ 5.99J/fc B , for D = 9 
we have T C (D) ~ 5.92J/&B- Since lateral phase separation 
occurs only for T < T C (D) [HHEH^], of course ' this shift 
restricts the range of £ that can be studied for the present 
choices of D. While 2£ ~ 1.0 occurs for T = 4.75J/&B and 
hence this choice, for which the cell size of the Ginzburg 
Landau model agrees with the lattice spacing of the molec- 
ular field model, is safely accessible, already for £ ~ 2.0 
(occurring for k^T — 5.875J) we are only slightly below 
T C (D = 9), and for £ = 3 we would already be in the one- 
phase region of such a thin film. In view of these consider- 
ations, k^T = 5.875J was chosen as a temperature where 
it makes sense to compare the lattice model with D = 9, 
19, and 29 (i.e., n m ax = 10, 20, and 30 lattice planes, re- 
spectively) with the corresponding Ginzburg Landau (GL) 
model. Note that for the convenience of comparison as well 
as accuracy of numerical solutions, the discrete mesh sizes 
Ax, Ay, and Az in the GL model were adjusted to the 
lattice constant a and we have set tq = 0. 

In Fig. [7] we present the comparison between the lat- 
tice model and GL model at k B T = 5.875J (£ ~ 2.0, 
m h ~ 0.25) with H S i = H S2 = 0.1 for D = 9 and 19. 
At this temperature the parameters hsi, 9, 1, T have the 
values 52.3 (523H S i), -128, 32, and t/1600, respectively. 
One sees that for D = 9 [Fig. E£a)] the behavior of the 
lattice model and the GL model are qualitatively similar, 
but there is no quantitative agreement. These discrepan- 
cies do get smaller, however, with increasing film thickness 
[Fig. EJb)]. In Fig- US where we plot the order parameter 
profile for D = 29 corresponding to two values of surface 
fields (see caption for details) , we see that the discrepancies 
between the lattice and GL models are already quite small, 
irrespective of the choice of the surface field. However, for 
stronger surface field [Fig. Hfb)] , the discrepancy is still vis- 
ible at the surfaces. We expect that the comparison should 
be perfect for any surface field in the semi-infinite limit. In 
the immediate vicinity of the critical point, the close cor- 
respondence between the time evolution predicted by the 
lattice theory and the GL model is also evident when one 
compares snapshot pictures of the time evolution, shown 
in Figs. HHTUl generated at k^T/J = 5.875 corresponding 
to Hsi = H S2 = 0.1, D = 29, for both the models. How- 
ever these snapshots suggest that at this temperature the 
lateral inhomogeneity disappears in the late stages of the 
phase separation process. Note that due to the proxim- 
ity of the critical point, which for finite D is shifted and 
nonzero surface fields, the coexistence curve separating the 
two-phase region from the one-phase region may be consid- 
erably distorted [l^,[3l|. In view of that, it is plausible that 
the state point for small D falls in the one-phase region of 
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FIG. 8: Same as Fig. [7] but for D = 29, and two choices of the 
strength of the surface fields, Hsi — Hs2 = 0.1 (a) and 1.0 (b). 



the thin film. However, this did not happen in the present 
case, as is clear from Fig. [11] where we have shown snap- 
shot pictures over a much longer time scale for a system 
with all parameters same as in Figs. HO [TU] except now we 
have D = 19, L — 64. So, the apparent stratified structure 
in Figs. [9j [10] is temporary which disappears at later time. 

In Fig. [T^] we show comparisons at lower temperatures, 
viz., k B T = 5.57J (C ^ 1.0, m b ~ 0.45, hsi ^ 15.4ff S i, 
g = -8, 7 = 4, t ~ i/104) and k B T ~ 4.75J (£ ~ 0.5, 
m b ~ 0.72, ~ 0.44Er S i, 5 = -0.5, 7 = 0.5, r ~ t/7.6). 
For both the temperatures we have set Hsi = = 0.1. 
Rather pronounced discrepancies between the lattice model 
and the GL model do occur, however, at low temperature 
[Fig. fl~2Tb)]. as expected. 

Note that the prefactor 7 [see Eq. ((2J] in the scaled 
Eq. {J} is an approximation which applies in the close vicin- 
ity of the critical point [8J . Here we try to take into account 
correction terms to the leading behavior of Eq. ([7]) to make 
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FIG. 9: Cross-sectional snapshot pictures in the a;z-plane, for 
D = 29, Hsi = H S 2 = 0.1, at k B T/J = 5.875 [same system as 
in Fig[8ja)], for four different times as indicated, according to 
(a) the lattice model and (b) the Ginzburg Landau model. 



the GL model more accurate for temperatures away from 
criticality. The order parameter <p(p,z,t) (not normalized 
by mb, and lengths not rescaled by 2£) satisfies the bound- 
ary conditions Q 



2r. 



d<t>{p, z = 0,t) 
dt 



T 

T c b 
T 



T \ J 

-1- J - 
T 



5j cj)(p,z = 0,t) 
d<j)(p,z,t) 



dz 



o,(29) 
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FIG. 10: Same as Fig. [9] but for a plane at n 
the walls. 
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15, parallel to FIG. 11: Same as Fig. [TOj but for D = 19, L = 64, at n = 10. 



+ ^^2^(P.*.*)]}|»=0 = 0. 



(30) 



From Eq. we see, however, that a pathological be- 
havior occurs if the coefficient (T c b/T — 1 — J/T) vanishes, 
which is the case for k^T/ J = 5: the time evolution of 
4>(p, z = 0, i) then is strictly decoupled from the order 
parameter in the interior, and it stops if 4>(p, z = 0, t) 
reaches the value Hsi/J (for J s /J = 1). This is what 
is seen in Fig. [TSTa') where we have solved the unsealed 
version of the GL model. In this case ^^(z = 0, t) has 
stopped its time evolution already during the very early 
stages. For k-gT < 5 J [Fig. [TBT b)]. the coefficient of 
the last term on the right side in Eq. has changed 
its sign (in comparison to the region close to T c b), and 



this leads to the result that 4>(p,z = 0, t) converges to 
zero, which also is unreasonable. Thus, taking the co- 
efficient (T ch /T - 1 - J/T) rather than simply {-J/T) 
[the latter leads to the scaled form Eq. ([7]) with the co- 
efficient 7 as quoted in Eq. does not yield any im- 
provement, but rather is physically inconsistent. However, 
working with the scaled form of the GL equations, and their 
boundary conditions, Eqs. ([U El El El does not yield re- 
sults in agreement with the lattice model at temperatures 
k B T/J < 5 either. 



V. CONCLUSION 

In this paper we have presented a Molecular Field The- 
ory for the Kawasaki spin-exchange Ising model in a thin 
film geometry and have shown that the numerical solution 



11 




FIG. 12: Plot of # a v(n) vs. n for D = 29 with Hsi = H S2 = 
0.1 at temperatures (a) k B T = 5. 57 J and (b) k B T = 4. 75 J. 
The open symbols corresponds to the lattice model whereas the 
filled symbols are for the GL model. Note that k B T = 5.57J 
corresponds to bulk correlation length £ ~ 1.0 and at k B T — 
4.75J, £ ~ 0.5. 



of the resulting set of coupled ordinary differential equa- 
tions describing the time-dependence of the local magne- 
tization at the lattice sites is a convenient and efficient 
method to study spinodal decomposition of such systems, 
taking the boundary conditions at the surfaces of the film 
properly into account. Obviously, in comparison to a 
Monte Carlo simulation of this model one has lost ther- 
mal statistical fluctuations, except for those built into the 
theory via the choice of noise in the initial configuration of 
the system; but there is consensus [H, HI H, 0, [B] that for a 
description of the late stage coarsening behavior such ther- 
mal fluctuations may safely be neglected. Thus, the present 
method is favorable in comparison with Monte Carlo, since 
the code runs much faster. 



FIG. 13: Comparison between the unsealed GL model [cf. 
Eq. (|29l30fl ] and lattice model at temperatures (a) 5J/k B and 
(b) AJ/k B . Note that the open symbols are for lattice model 
and the filled symbols are used for the unsealed GL model. 



In the vicinity of the critical point, where the correla- 
tion length £ is sufficiently large, and also the film thick- 
ness D is sufficiently large as well, our treatment becomes 
equivalent to the time-dependent Ginzburg-Landau theory. 
However, one needs to go surprisingly close to the critical 
point of the bulk to actually demonstrate this limiting be- 
havior from our lattice model treatment numerically. The 
GL treatment over most of the parameter regime provides 
only a qualitative, rather than quantitative, description of 
the system. In principle, for the GL theory to be accurate, 
the correlation length should be much larger than the lat- 
tice spacing. However, this happens only very close to the 
critical point. Away from the critical points no well-defined 
connection to the parameters of a microscopic Hamiltonian 
can be made for the GL theory, while the present lattice 
approach has this connection by construction. 
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Of course, the problem of surface-directed spinodal de- 
composition is most interesting for liquid binary mixtures, 
for which our lattice model is inappropriate due to the 
lack of hydrodynamic interactions; even for solid mixtures 
(two different atomic species sharing the sites of a lattice) 
our model is an idealization (neither elastic distortions nor 
lattice defects were included; actual solid binary mixtures 
decompose via the vacancy mechanism of diffusion [l|, Q ; 
etc.). While GL models including hydrodynamic interac- 
tions have been formulated p], [H, d, d, [f| , the restriction 
that any GL model is valid only in the immediate neigh- 
borhood of the critical point applies there as well. While 
outside the critical region unmixing of fluids can be simu- 



lated by Molecular Dynamics methods (see e.g. [Hj]), such 
simulations are extremely time-consuming, and an analog 
of the present dynamic mean field theory for inhomoge- 
neous fluids would be very desirable. Developing such an 
approach clearly is a challenge for the future. 
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